#pragma once

#include "utility.h"
#include "math.h"
#include <vector>



//#include "Vector.h"

//+==================================================================+
//////////////////////////////////////////////////////////////////////
//	dSFMT extension fors SSE datatypes
//+------------------------------------------------------------------+
inline __m128d dsfmt_genrand_close1_open2_pd(dsfmt_t *dsfmt)
{
	double *psfmt64 = &dsfmt->status[0].d[0];

// make index a multiple of 2, for alignment reasons
	int &i = dsfmt->idx;
	i = (i+1) & ~1;

	if(i+1 >= DSFMT_N64) {
		dsfmt_gen_rand_all(dsfmt);
		i = 0;
	}

	__m128d r = *(__m128d*)(psfmt64+i);
	i += 2;
	return r;
}


//+------------------------------------------------------------------+
inline __m128 dsfmt_genrand_close1_open2_ps(dsfmt_t *dsfmt)
{
	double *psfmt64 = &dsfmt->status[0].d[0];

// make index a multiple of 2, for alignment reasons
	int &i = dsfmt->idx;
	i = (i+1) & ~1;

	if(i+3 >= DSFMT_N64) {
		dsfmt_gen_rand_all(dsfmt);
		i = 0;
	}

	__m128d *src = (__m128d*)(psfmt64+i);
	__m128 r = _mm_movelh_ps(_mm_cvtpd_ps(src[0]),_mm_cvtpd_ps(src[1]));

	i += 4;
	return r;
}

//+------------------------------------------------------------------+
inline void dsfmt_fill_array_close1_open2_ps(dsfmt_t *dsfmt, __m128 *dst, unsigned count)
{
	double *psfmt64 = &dsfmt->status[0].d[0];

	__m128 *lim = dst + count;

	unsigned step = DSFMT_N64;
	step >>= 2;

	for(;;)
	{
		dsfmt_gen_rand_all(dsfmt);

		if(dst + step <= lim)
		{
			cvtpd2ps(psfmt64, (float*)dst, step*4);
			dst += step;
		}
		else
			break;
	}

	ptrdiff_t remainder = lim - dst;
	if(remainder)
	{
		remainder *= 4;
		cvtpd2ps(psfmt64, (float*)dst, remainder*4);
		dsfmt->idx += remainder;
	}

}


//+==================================================================+
//////////////////////////////////////////////////////////////////////
class RNG
{
	dsfmt_t dsfmt;

public:
//+------------------------------------------------------------------+
	RNG(unsigned int s = 0)
	{
		if(s == 0)
		{
			union {u64 a; u32 b[2];} data;
			data.a = __rdtsc();
			dsfmt_init_by_array(&dsfmt,data.b,2);
		}
		else
			seed(s);
	}

//+------------------------------------------------------------------+
	void seed(unsigned int s)
	{
		dsfmt_init_gen_rand(&dsfmt,s);
	}

//+------------------------------------------------------------------+
	static int GetMinArraySize()
	{
		static int minArraySize = dsfmt_get_min_array_size();
		return    (minArraySize);
	}

//+------------------------------------------------------------------+
//	A uniformly distributed random number from the range [1, 2)
	double rand12()
	{
		return dsfmt_genrand_close1_open2(&dsfmt);
	}

//+------------------------------------------------------------------+
//	A uniformly distributed random number from the range (0, 1)
	double rand_open()
	{
		return dsfmt_genrand_open_open(&dsfmt);
	}

//+------------------------------------------------------------------+
//	A uniformly distributed random number from the range [0, 1)
	double rand()
	{
		return dsfmt_genrand_close_open(&dsfmt);
	}

//+------------------------------------------------------------------+
//	A uniformly distributed random number from the range [0, 1)
	float randf()
	{
		return (float)rand();
	}

//+------------------------------------------------------------------+
//	Two uniformly distributed random numbers from the range [1, 2)
	__m128d rand12_pd()
	{
		return dsfmt_genrand_close1_open2_pd(&dsfmt);
	}

//+------------------------------------------------------------------+
//	Two uniformly distributed random numbers from the range (0, 1)
	__m128d rand_open_pd()
	{
		return _mm_castsi128_pd(_mm_or_si128(_mm_castpd_si128(rand12_pd()),sse2_int_one)) - _pd_1;
	}

//+------------------------------------------------------------------+
	__m128d rand_pd()
	{
		return rand12_pd() - _pd_1;
	}

//+------------------------------------------------------------------+
	__m128 rand12_ps()
	{
		return dsfmt_genrand_close1_open2_ps(&dsfmt);
	}

//+------------------------------------------------------------------+
//	Four uniformly distributed random numbers from the range [0, 1)
	__m128 rand_ps()
	{
		return rand12_ps() - _ps_1;
	}

//+------------------------------------------------------------------+
//	Uniformly distributed random numbers from the range (0, 1)
	__m128 rand_open_ps()
	{
		return _mm_castsi128_ps(_mm_or_si128(_mm_castps_si128(rand12_ps()),_epi32_1)) - _ps_1;
	}

//+------------------------------------------------------------------+
//	Uniformly distributed random integer from the range [0, limit)
	template<class T>
	T randInt(T limit)
	{
		return (T)(limit*dsfmt_genrand_close_open(&dsfmt));
	}

//+------------------------------------------------------------------+
//	0 or 1 with equal probability
	bool randBool()
	{
		return (dsfmt_genrand_uint32(&dsfmt) & 0x80000000U) == 0;
	}

//+------------------------------------------------------------------+
//	Generates any of the possible 32 bit integers with equal probability.
	unsigned int randomBits()
	{
//		return 4294967296.0 * rand();
		return dsfmt_genrand_uint32(&dsfmt);
	}

//+------------------------------------------------------------------+
	__m128 random_bits_ps()
	{
		M128 r;
		r.u32[0] = dsfmt_genrand_uint32(&dsfmt);
		r.u32[1] = dsfmt_genrand_uint32(&dsfmt);
		r.u32[2] = dsfmt_genrand_uint32(&dsfmt);
		r.u32[3] = dsfmt_genrand_uint32(&dsfmt);
		return r;
	}

//+------------------------------------------------------------------+
//	Array of uniformly distributed random numbers from the range [1,2)
	void randFillArray12(double arr[], int count)
	{
		dsfmt_fill_array_close1_open2(&dsfmt, arr, count);
	}
//+------------------------------------------------------------------+
//	Array of uniformly distributed random numbers from the range [0,1)
	void randFillArray(double arr[], int count)
	{
		dsfmt_fill_array_close_open(&dsfmt, arr, count);
	}
//+------------------------------------------------------------------+
//	Array of uniformly distributed random numbers from the range [1,2)
	void randFillArray12_ps(__m128 arr[], unsigned count)
	{
		dsfmt_fill_array_close1_open2_ps(&dsfmt, arr, count);
	}
//+------------------------------------------------------------------+
//	Array of uniformly distributed random numbers from the range [0,1)
	void randFillArray_ps(__m128 arr[], unsigned count)
	{
		randFillArray12_ps(arr, count);
		for(unsigned i = 0; i < count; i++)
			arr[i] -= _ps_1;
	}

//+------------------------------------------------------------------+
//	A gaussian distributed random number with mean=0 and std.dev.sigma
//	~10 times slower than rand().
	double gauss(double sigma)				{	return gauss() * sigma;	}
//+------------------------------------------------------------------+
	double gauss()
	{
		M128 g;
		g.m128d = gauss_pd();
		return g.f64[0];
	}

//+------------------------------------------------------------------+
//	Two gaussian distributed random numbers,
//	with zero mean and one std.dev.
	__m128d gauss_pd()
	{
		M128 x, w;

		do
		{
			x.m128d	= dsfmt_genrand_close1_open2_pd(&dsfmt);
			x.m128d *= _pd_2;
			x.m128d -= _pd_3;

			w.m128d	= _mm_mul_pd(x.m128d,x.m128d);
			w.f64[0]+= w.f64[1];
		}
		while(w.f64[0] >= 1.0);

		w.f64[0] = sqrt( -2.0*log(w.f64[0]) / w.f64[0] );
		w.f64[1] = w.f64[0];

		return _mm_mul_pd(x.m128d,w.m128d);
	}


//+------------------------------------------------------------------+
//	A cauchy distributed random number.
	double cauchy(double gamma = 1.0)
	{
		return gamma * tan(PI*(rand_open()-0.5));
	}

//+------------------------------------------------------------------+
//	An exponentially distributed random number.
	double exponential(double mu = 1.0)
	{
		return -mu * log(2.0-rand12());
	}

//+------------------------------------------------------------------+
//	Four exponentially distributed random numbers. Multiply by mu.
	__m128 exponential_ps()
	{
		return -log_ps(_ps_2-rand12_ps());
	}

//+------------------------------------------------------------------+
//	A levy distributed random number, with scale c and exponent alpha.
//	0 < alpha <= 2
//	437 cycles on my Core2Duo.
	double levy(double c, double alpha);

//+------------------------------------------------------------------+
//	Levy distributed random numbers, with scale c and exponent alpha.
//	0 < alpha <= 2
//	625 cycles on my Core2Duo.
	__m128 levy(__m128 c, __m128 alpha);

//+------------------------------------------------------------------+
//	Uses levy distribution to generate a random number between (lo,hi)
//	but usually close to the original number x.
	__forceinline double levyMutation(double x, double lo = 0.0, double hi = 1.0, double c = 0.35, double alpha = 0.35)
	{
		double xNew;
		do
		{	xNew = x + levy(c,alpha);
		}
		while(xNew >= hi || xNew <= lo);

		return xNew;
	}
//+------------------------------------------------------------------+
	__forceinline float levyMutation(float x, float lo = 0.0f, float hi = 1.0f, double c = 0.35, double alpha = 0.35)
	{
		float xNew;
		do
		{	xNew = x + (float)levy(c,alpha);
		}
		while(xNew >= hi || xNew <= lo);

		return xNew;
	}

//+------------------------------------------------------------------+
//	Fisher-Yates shuffle
//	n = the number of items left to shuffle
	template<class T>
	inline void shuffle(T *v, int n)
	{
		while(n > 1)
		{
			int k = randInt(n);					// 0 <= k < n.
			n--;										// n is now the last pertinent index;
			std::swap(v[n],v[k]);				// swap array[n] with array[k] (does nothing if k == n).
		}
	}

//+------------------------------------------------------------------+
	template<class T>
	inline void shuffle(std::vector<T> &v, int n = -1)
	{
		if(n < 0)
			n = v.size();

		return shuffle(&v[0],v.size());
	}

};

